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CN ■ Abstract 



The Van Leer approach for the approximation of nonhnear scalar conservation 
[ laws is studied in one space dimension. The problem can be reduced to a nonlinear 
interpolation and we propose a convexity property for the interpolated values. We 
prove that under general hypotheses the method of lines in well posed in n BV 
and we give precise sufficient conditions to establish that the total variation is 
diminishing. We observe that the second order accuracy can be maintained even 
at non sonic extrema. We establish also that both the TVD property and second 
order accuracy can be maintained after discretization in time with the second 
^ ■ order accurate Heun scheme. Numerical illustration for the advection equation is 
m ; presented. 
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Interpolation non lineaire 
et schemas a variation totale decroissante 



Nous etudions le schema de Van Leer pour I'approximation de lois de conservation 
scalaires a une dimension d'espace. Nous proposons une propriete de convexite 
de I'interpolation non lineaire associee a cette methode. Nous prouvons que sous 
des hypotheses generales, le schema continu en temps conduit a un probleme bien 
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pose et a variation totale decroissante. Ce resultat s'etend au cas discretise en 
temps et des tests numeriques mettant en evidence divers ordres de convergence 
sont presentes. 
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1. Introduction 

We consider the following scalar conservation law : 

Ou 

(1.1) — + — /(w) = , xeJR, t>0, u(x,t)eJR 
at ox 

where the flux function / is supposed to be a convex regular (C^ class) real 
function. The initial condition takes the form 

(1.2) u{x, 0) = uo{x) , xeU 

and we assume that the initial datum uq satisfies : 

(1.3) uo€L°°(IR) n BV(IR). 

• Even for smooth initial datum wq, a classical solution of the Cauchy problem 
(1.1). (1.2) does not exist in general {e.g. SmoUer [1983]) and we consider in this 
paper weak solutions of the problem (1.1). (1.2), i.e. functions u G L°°(IR x 
[0, +oo[) such that 
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(1.4) 




V^gCo(IRx]0, +oo[) 



and 



(1.5) 



t 



lim 
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• The uniqueness of a weak solution satisfying (1.4). (1.5) is guaranteed if 
we consider (Lax [1971]) only entropy solutions of the conservation law (1.1), i.e. 
weak solutions also satisfying 

(1.6) [ \r]{u)^ +au)^]dxdt <0, G Co^KxlO, +oo[) , ^ > 



dt dx 



where r] : JR, i — > IR is a strictly convex entropy and ^ : IR i — > IR is the 
associated entropy flux : 

(1.7) e{u) = V'{u).f'{u). 
With the particular choice 

(1.8) Vc{u) =\u — c\, ceIR, 

Kruskov [1970] proved that a weak solution of (1.1). (1.2) satisfying the inequalities 
(1.6) for all the choices of entropies (1.8) is necessarily unique. More recently, Di 
Perna [1983] proved that uniqueness is implied if the inequality (1.6) is valid for a 
single strictly convex entropy r]. 

• The approximation of (1.4)- (1.6) by finite volume numerical schemes is a 
very classical problem. In this paper, we first look at approximations that are 
discrete in space and continuous in time (method of lines) . Let h > be a real 
parameter and xj = j h he the associated mesh points. Following Godunov [1959], 
we search an approximation Uj{t) which is constant in the interval 



1 rU+h)h 



(1.9) Uj{t) « T / u{x,t)dx, jelZ, t>0. 

h J{j-^)h 

As was recognized by Lax-Wendroff [1960], we require Uj{t) to satisfy a conser- 
vative numerical scheme : 

^^■^^^ IF + 1(^^+4 -■'■^-0 = °' ^'=^' 

• The system of ordinary differential equations (1.10) is now completely de- 
fined if we prescribe the initial conditions and the numerical flux /j+i. Three 
point numerical schemes are parameterized by the numerical flux function $ : 

(1.11) /j+i = ^{Uj,Uj+i) 

which is supposed to be Lipschitz continuous and consistent with the flux / : 

(1.12) ^{u,u) = f{u), w € IR. 

The numerical flux realizes some approximation of the Godunov flux (we refer 
to Harten-Lax-Van Leer [1983] for the details). Following Crandall-Majda [1980], 
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the flux / of relation (1.11) is monotone if (w, v) i — > v) is a nondecreasing 
(respectively nonincreasing) function of u (resp. v). Kusnezov-Volosin [1976] and 
Crandall-Majda [1980] have extended to monotone schemes the proof of conver- 
gence towards the entropy solution developed by Leroux [1976] for the Godunov 
scheme (when the equation (1.10) is also discretized in time). Sanders [1983] 
proved that on irregular meshes, the solutions of the ODE (1.10) also converges 
to the entropy solution of (1.4)-(1.6) as h tends to zero. Osher [1984] with the 
"E-schemes" extended the notion of monotone schemes to take into account non 
necessarily convex fluxes. The major default of three point monotone schemes is 
that they can have at best flrst order accuracy. 

• High resolution schemes (that are in general flve point numerical schemes) 
have been proposed to overcome this lack of precision. On one hand, Harten 
[1983] introduced the notion of total variation diminishing (TVD) schemes ; if the 
discrete total variation 

(1.13) TV(t) = Yl \^j+iit)-Mt) I 

is a nonincreasing function of time, then the scheme (1.10) is said to be TVD. 
Monotone schemes are particular cases of TVD schemes and Harten gave sufficient 
conditions (generalized by Sanders [1983] for the method of lines) to ensure that 
a numerical scheme is TVD. Unfortunately, three point TVD schemes are at most 
first order accurate and Harten [1983] constructed second order numerical schemes 
by modifying the first order fiux (1.11) into a five point scheme in a way analogous 
to the Flux Corrected Transport (FCT) of Boris-Book [1973]. 

• On the other hand, the "Multidimensional Upstream centered Scheme for 
Conservation Laws" (MUSCL) approach proposed by Van Leer [1979] (see also 
Colella [1985]) assumes that a piecewise linear interpolation is reconstructed from 
the mean values Uj (t) : 

(1.14) Uh{x, t) = Uj{t) + Sj{t) {x — Xj) , X Elj , j E . 

The slopes Sj are chosen as nonlinear functions of the three mean values at the 
points j - 1, j, j + 1 : 

(1.15) Sj = S{uj-i , Uj , Uj+i) 

in such a way that the reconstruction (1.14) satisfies some monotonicity restrictions 
(Van Leer [1977], see also Sweby [1984] for a unified vision of the fiux limiters). 
The function x i — y Uh{x, t) is a priori discontinuous at the points x^,i ; the 

two values u^, i are defined on each side of x,, i : 
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(1.16) (a) u._^, = Uj{t) + -S{uj-i, Uj , Uj+i) 

(1.16) (b) u+^i = Uj+x{t) - ^ S{uj , Uj+i , Uj+2) ■ 

The flux function /j+i is computed as in the Godunov scheme, but with the 
latter two values as arguments : 

(1.17) /i+i = ^^0' 

The theoretical convergence properties of MUSCL schemes have been studied by 
Osher [1985]. He gave sufficient conditions to establish that the scheme (1.10). (1.17) 
is TVD and proved under some restrictions the convergence to the entropy solution 
(using a discrete entropy inequality proposed in Osher [1984]). Moreover second 
order accuracy is realized even at regular extrema which are not sonic points, i.e. 
values of u such that 

(1.18) f'{u) ^ 0. 

This fact was recognized by Osher- Chakravarthy [1984]. The MUSCL approach 
has been extended to higher orders by Collela- Woodward [1984] (PPM method) 
and by Harten et al [1987] with the so-called ENO schemes. 

• The main purpose of this paper is to show how to choose nonlinear inter- 
polation formulae like (1.16) in order to get a TVD scheme. In the second part of 
this paper we propose some natural limitations on the interpolated values between 
two mesh points and we establish in the third part sufficient conditions that prove 
that the resulting scheme is TVD. In part IV, we demonstrate that the second 
order accuracy is obtained in all smooth parts of the flow. In part V we dicretize 
the method of lines in time and we precise conditions to maintain both the TVD 
property and second order accuracy in smooth regions. In the last part we present 
some numerical experiments for the advection equation with periodic boundary 
conditions. 



2. Some Limitations for the Interface Values 

We recall that a sequence of values Uj is supposed to be given at the points Xj of 
a regular mesh on the real axis. We notice a small difference with the finite volume 
method where the variables Uj represent the mean values of some function in each 
interval Ij. We interpolate those values at the points x^j^i by two different values 

u~, 1 and M^, 1 on each side of the interface , 1 . By analogy with the MUSCL 

approach (1.16), we set 
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(2.1) u . = L{uj-i, Uj , Uj+i) 

2 

(2.2) w+ 1 = R{uj , Uj+i , Uj+2) ■ 

2 

The interpolation defined by (2.1). (2. 2) is a priori nonlinear and in the follow- 
ing of this part we give some restrictions concerning the functions !/(•, •, •) and 



Property 2.1 Homogeneity 

We suppose that the multiplication of all the arguments of (2.1). (2. 2) by the same 



ment, i.e. that the functions L and R are homogeneous of degree 1 : 

(2.3) J{Xu,Xv,Xw) = XJ{u.,v,w), J = LovR, AelR. 

Property 2.2 Translation invariance 

We suppose that by addition of some constant A to the three arguments of 
(2.1). (2. 2), the resulting interpolate values are translated by the same value : 

(2.4) J{u-\- X, V -\- X, w -\- X) = J{u,v,w) + X, J = LoxR, AelR. 
Property 2.3 Left-Right symmetry 

We suppose that by exchange of the two extreme arguments of (2.1) {i.e. we 
reverse left and right) the left interpolation at x,,i is transformed into the right 
interpolation at Xj_i (Figure 2.1) : 

(2.5) R{w., v., u) = L{u, V, w) . 



R{», •, •). 



real number A multiplies the interpolate values u . 



3 + ^ 



w 1 by the same argu- 



u — 




u 



j-1 



J 



j+1 



j-1 



J 



j+1 



Figure 2.1. Left-right symmetry. 
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□ 



Proposition 2.1 

If the properties 2.1 to 2.3 are assumed, there exists some hmiter function ip 
JR — y JR such that 

(2.6) L[u, V, w) = V -\- -yjl ] [v — u) 

2 \v — u J 

1 / v — u\ 

(2.7) R{u^ V, w) ~ ~ 2 ^\ / i'^ ~ ■ 

Proof of Proposition 2.1 

We have simply 

/ w — V \ 

Liu, V, w) = V -\- Liu — V, 0, w — v) = V -\- (v — u) Li — 1, 0, ) 

V V — u J 

and R{u^ v, w) = L(w, u) = + 2 ^( ) ("^ ~ '"^) • 

Property 2.4 Monotonicity 

We suppose that the interpolated values u~_^i and u^_^i belong to the interval 

[min(wj, Wj+i) , max(wj, Wj+i)], i.e. 

(2.8) min(i;, w) < L{u, v, w) < max('y, w) 

(2.9) min(w, v) < R{u, w) < max(w, v) . 

We remark that the inequalities (2.9) are a simple consequence of (2.5) and (2.8). 
Proposition 2.2 

If the properties 2.1 to 2.4 are valid, the limiter ip satisfies the condition 

(2.10) < ^ < 2, AgIR 

A 

and if the properties 2.1 to 2.3 are supposed, the inequalities (2.10) imply the 
property 2.4. 

Proof of Proposition 2.2 

It is a straightforward consequence of (2.6), written under the form : 

iH—) 

(2.11) L{u,v,w) = v + - (w-v). □ 

V — u 

• We introduce now an original notion concerning the nonlinear interpolation 
at the interfaces, first proposed in Dubois [1988]. 



FRANgois Dubois 



Property 2.5 Convexity 

We suppose that for each integer j, the sequence of the interpolate values 

(2.12) Uj-i , w+i , Uj , w7 1 , Uj+i , 

■J 2 J' 2 

is the restriction to the nodes 

(2.13) Xj-i 5 -^j'-i ' -^j ' ' "^j+i ' 

of some convex or concave function, depending on the concavity of the sequence 

(2.14) Uj_i , Uj, , Uj+i . 
More precisely, if the inequality 

(2.15) Uj — Uj-i < Wj+i — Uj 

holds (convexity of the sequence (2.14)), we suppose that the discrete increments 
of (2.12) realizes a nondecreasing sequence : 

(2.16) w+i - Uj-i < Uj - w+i < w7 1 - Uj < Uj+i - w~ 1 

J 2 J 2 J' 2 J' 2 

and if we change "<" into ">" in the relations (2.15) then we must do it also in 
(2.16) (concavity of the sequence (2.16)), see Figure 2.2. 



A 



1 \ 
1 \ 
1 

1 




j-1 


j-1/2 


j 


j+1/2 


j+1 



Figure 2.2. Convexity of the interpolated value at the intermediate points. 
Example 2.1. 

The classical Lagrange interpolation by a polynomial of degree 2 inside the cell 
]Xj_i, Xjj^i[ corresponds to the limiter function = ^ . It satisfies the 

properties of homogeneity, translation invariance, convexity but not monotonicity 
(see Van Leer [1973]). 
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Property 2.3 

If the properties 2.1 to 2.3 hold then the property 2.5 is equivalent to the following 
inequalities concerning the limiter function : 

(2.17) 1 < A^(i) < ^(A) < A, A>1 

(2.18) A < V^(A) < AV^(i), A<0. 
Proof of Proposition 2.3 

We first suppose that the property 2.5 is valid. We set 

(2.19) w = V + X{v-u). 

Then the inequality (2.15) is satisfied if and only if 

(2.20) (A-1) (v-u) > 0. 

Then the inequalities (2.16) take the form 

(l-^X^J;iX)){v-u) < -^X^I;{X){v-u) < ... 

... < ^i^{X){v-u) < (a - ^(A)) (i; - n) . 

If we suppose that [v — u) > 0, we obtain easily (2.17). If we reverse the inequality 
(2.20), we have also to reverse (2.21), thus (2.17). If < A < 1 the change of 
variable M = ^ gives 

1 > - V^(/x) > ^(-) >-, A=i, 0<A<1 

which is exactly (2.17). When A < 0, the inequality (2.21) remains valid but 
under the hypothesis v — u < 0. Then we have 

(2.22) 1 > AV^Q) > V^(A) > A, A<0 

which establishes (2.18). 

• We suppose now that the inequalities (2. 17). (2. 18) hold. We claim that 
(2.22) is valid because if we change A into j in the inequality (2.18) we get 



(2.21) 



i < < i^(A), A<o 



which establishes the two first inequalities of (2.22). The end of the proof can be 
obtained without difficulty. □ 



We summarize the last two properties (relations (2.10), (2.17), (2.18)) into the 
following one : 
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Proposition 2.4 

If we suppose that the properties 2.1 to 
(2.6) satisfies 

(2.23) 1 < A ^(i) < ^(A) < A 

(2.24) A < ^(A) < 2A 

(2.25) A < ^(A) < 



2.5 are valid, the limiter t/j defined in 

A > 1 

< A < 1 
A < 0. 



Figure 2.3 represents the constraints associated with (2. 23). (2. 25) ( except the 
second inequality of (2.23) !). 




Figure 2.3. Constraints of monotonicity and convexity for the limiter function 
t/j defined at (2.6). The graph of A i — > tpiX) lies inside the shaded region. 



Remark 2.1. 

The interpolated values u~_i and satisfy (according to (2.1), (2.2), (2.6), 

2 2 

(2.7)) : 

(2.26)(a) = uj + - i/j{Xj) [uj — uj-i) 



(2.26) (b) 
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and following e.g. Sweby [1984] we can introduce the inverse rj of the ratio Xj 



(2.27) Tj = 

Uj + i - Uj 

and rewrite the limiter ip in terms of the variable r, by defining 



(2.28) (/p(r) = r^(i) 



We emphasize that there is a priori no reason to suppose = i/j. The two first 
equalities of (2.26) take the form 

(2.29)(a) uT^, = uj + ^<f{rj){uj+i-Uj) 

1 1 

(2.29) (b) u+_, = Uj - -^(-) (uj-Uj^,). 

Moreover, the inequalities derived in Proposition 2.4 can be translated in terms of 
the limiter cp defined in (2.28). 

Proposition 2.5 

Let (p be defined by (2.28). The limiter i/j satisfies the inequalities (2.23)- (2.25) 
if and only if the function c/? verifies : 

(2.30) max(0, r) < ip{r) < 1, -oo < r < 1 

(2.31) 1 < ip{r) < min(r, 2), r > 2 

(2.32) (p{r) < r (p(^) , r > 1. 
Proof of Proposition 2.5 

• We first suppose r < 0. Thus the inequality (2.25) joined with (2.28) shows 
that we have < (p{r) < 1 and (2.30) is established in this case. If r > 0, the 
convexity property reduces to (2.17) and the monotonicity property corresponds 
to (2.10). By changing A into ^, we get 

(2.33) r < r <f{^) < (^(r) < 1, < r < 1 

which achieves the proof of (2.30). 

• If r > 0, we get from (2.33) : 

(2.34) 1 < ip{r) < r, < r > 1 

and also (2.32). On the other hand condition (2.10) can be rewritten by changing 
A into - ; we deduce : 

(p(r) = r ip(^-^ < 2, r > 1. 
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The latter inequality, joined with (2.34) proves that (2.31) holds. The other side 
of the equivalence is established by similar arguments. □ 

Example 2.2. 

Let ijj be defined by the relations 

(2.35)(a) ^|J{X) = 0, A < 

(2.35)(b) ^(A) = 2A, 0<A<i 

(2.35)(c) ^{X) = ^ + ^A, ^<A<3 

(2.35) (d) ^|J{X) = 2 , A > 3. 

This corresponds to the original MUSCL extrapolation (Van Leer [1977]) and the 
K limiters af Anderson, Thomas, Van Leer [1986] (see also Chakravarthy et al 
[1985]) correspond to the choice 

(2.36) (a) V^,(A) = 0, A < 

(2.36) (b) ^,(A) = i I (1-/.) min (l , ^ a) + (1+/.) min (a , | 

if A > 0. All of those limiters (—1 < k < 1) satisfy ip = if) (see (2.28)) and 
the inequalities (2.23). (2.25) ; we remark that in the latter relations, we do not 
impose the values of ip corresponding to negative arguments to be null. This last 
condition is necessary if we interpolate u by a linear function inside the interval 
Ij and suppose moreover the condition of monotonicity. We remark also that the 
so-called Van Albada limiter [1982], defined by 

(2.37) ^(A) = :|^1±^ = ^^G)' 

satisfies neither (2. 23). (2. 25) [for A < -1] nor (2.30). (2.32) [for -1 < A < 0]. 
Proposition 2.6 

The limiter L{u,v,w) satisfying the properties 2.1 to 2.3 is Lipschitz continuous, 
i.e. 

3C > 0, Vw, w, u\ v\ w' e IR, 

I I/(w, V, w) — L{u', v', w') |< C (y \ u — u \ + \v — v' \ + \ w — w' \ ) 

if and only if the limiter functions ip (defined at (2.6)) and (defined at (2.28)) 
satisfy a similar property : 

(2.39)(a) 3i^>0, VA, A'gIR, | V^(A') - V^(A) | < | A' - A 

(2.39)(b) 3X > 0, Vr, / € IR, | (^(r') - <^(r) | < K\r'-r\ 



(2.38) 
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Proof of Proposition 2.6 

• We suppose first that (2.38) holds. We take u = —1, v = 0, w = X and 
analogous relations for the "prime" variables. Then we have (according to (2.6)) : 

i (^{X') - V^(A)) = L(-l, 0, X'-l) - L(-l, 0, A - 1) . 
Similarly, we have : 

i {ip{r') - ^(r)) = L(-/, 0, 1) - L(-r, 0, 1) 
which establishes (2.39). 

• We suppose that (2.39) is realized. For a set of given real numbers u, w 
(respectively u', v', w') , we define A, r (resp. A', r') by the formulae w = v + 
X{v — u) ] V = u-\- r{w — v) (resp. w' — v' + X'{v' — u') \ v' = u' + r'{w' — v')) and 
we decompose the left hand side of (2.38) as follows : 

( L{u' , v', w') — L{u, V, w) = v' — v-\- 
^^■^^^ I + ^ (^(rO - ^(r)) {w' -w) + ^ ip{r) {{w' - w) - {v' - v)) 

If r' = oo (resp. r = oo) then w' = v' (resp. w = v) and (2.38) follows directly 
from (2.40) and (2. 30). (2. 31). In the other cases, the estimate 

(2.41) \r'-r\\w'-v'\<\u'-u \ + (1+ | r |) \v' -v \ + | r | | k;' - k; | 
joined to (2.39) establishes the property (2.38) under the restriction 

(2.42) I r I < 1. 

The same conclusion holds if we reverse the roles of r and r' . Inequality (2.38) is 
established if we suppose 

(2.43) I A I < 1 and | A' | < 1 . 

We have in this case 

{L(m', v' , w') — L{u, V, w) — v' — v-\- 
+ i ^(A) {v'-v + u'-u) + ^ (V^(A0 - ^(A)) {v' - u') 

and according to (2.10) and (2.39) it is sufficient to remark that we have 

I A' - A I \v' -v \ < (1+ I A I) \v' - v\ + I A I I - u I + | w;' - w; | 
to obtain finally : 

I L(w', v', w') - L{u, v,w)\< {l+'^^{\u' -u \ +2 \v' -v \ + \ w' -w\) 
which ends the proof. □ 
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+ 
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Remark 2.2. 

The total variation of the "reconstructed sequence" 

(2.45) u = ■■■ , Uj-i , w7_ 1 , w+_ 1 , Uj , w" 1 , w+ i , w^+i , • • • (j e ^) 

J 2, J 2 J' 2 J' 2 

is naturahy defined by 

(2.46) TV(w) = J] (I I + I ^ 

The properties 2.1 to 2.5 do not imply that the total variation of the reconstructed 
sequence is smaller than TV(u) (c./. (1.23)). For example, in the particular case 
of the following simple sequence 

(2.47) (a) Wj- = -3/i, J < -1 
(2.47)(b) Wo = 0, ui = h 

(2.47) (c) Uj = 4h, j >2 

the reconstructed sequence associated with the MUSCL scheme (2.35) admits a 
total variation given by 

(2.48) TV(w) = TV(w) + 2/i. 
(see Figure 2.4) 



□ • □ • □ 




□ • □ • □ 



Figure 2.4. The total variation of the reconstructed sequence (2.38) with the 
MUSCL scheme (squares) is greater than the TV of the initial sequence (points). 
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3. Decrease of the Total Variation 

In this section, we suppose that the interpolated values at the interfaces 
(i + i)^ given by the relations (2. 27). (2. 29) and are submitted to the re- 
strictions described in the Part II (properties 2.1 to 2.5). We study the following 
differential system 

IT + 1 - ^^-i) = *>o 

with the flux at the interface {j + ^) h given by a monotone flux function $ : 

(3.2) v., «;^.), jez. 

We suppose that the initial datum 



(3.3) 


uj{0) = 


satisfies 




(3.4) 




id est 




(3.5) (a) 


sup 






(3.5) (b) 


El 







u° I < OO 



Uj+i — \ < oo 



Then the method of lines is well defined, at least locally in time. 
Proposition 3.1 

The semi-discrete numerical scheme defined by the differential system (3.1). (3. 2) 
and (2. 27). (2. 29) associated with the initial condition (3. 3). (3. 4) has a unique 
solution Uj{t), j e !Z, t G [0,r], for some T > 0, if the fiux function <l> is 
Lipschitz continuous on the compact sets of IR : 

u\< M, \v \< M, \ u' |< M, 



r MM > 0, 3C > 0, Vw, u' , v' eWi, ( 
^ ^ \ \v' \<M) =^\^{u\v')-^{u,v)\ 



< C{\u-u' \ + \ v-v' \ ) 



Proof of Proposition 3.1 

• From the classical Cauchy-Lipschitz theorem it is sufficient to prove that 
the difference of the flux 

(3.8) F,{u) = -i i - f^_.) , i e ^, n G E^{7L) 

is a Lipshitz continuous i.°^{7L) — )■ £°°{7L) function. According to (3.5), we 
suppose that \uj |< M. Then the same property is realized for the intermediate 
states li^, 1 according to the property 2.4 (monotonicity) . We have : 
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Fj{u) I < - I ^{u.^, , J - , Uj) I + I , uj) - ^Uj , uj 



3 + 1 



J 2 J 1 J 1 J 1 



< 
< 



ll 

C r 



u 



+ 



+ u 



Un 



+ I w+ 



Uj \ -\- \ u 



Un 



1 / 1 ^ , 



1 , , , , 1.1. 

+ ^'^{rj) I Uj - Uj-i I + I 1 - 2'^(^) 



Uj 



3-2 
Uj 

■ u 



Un 



C 

< 16 — M 

h 

that establishes that F is a function £°°{2L) — y i°°{2L). 

• We prove that F is Lipschitz continuous on the bounded sets of £°°{7L). 
Let {uj)jQ2Z: ivj)j^2L be two bounded sequences in £°°(7L) : 

(3.9) \uj\< M, \vj\< M, Vj€^. 

According to the monotonicity property (2.8), the interpolated values w^, i , 
I , satisfy the same property. From the inequality (3.7) and the proposition 

3~T~ 2 

2.6 we deduce 

I . 4+ J - i>(v. . ) I < ( I vi - vi I + 1 - I ) 

< KC (^\ Uj-i — Vj-i I +2 \uj—Vj\ +2 I Wj+i— i'j+i I + I Uj-\-i — Vj-\.i \ ) 

< QK C sup I — i;^ I . 

Then we have the estimate 

12 C 

(3.10) I Fj{u) - Fj{v) I < — - — sup \ui-Vi\ 
and the proposition is established. 



□ 



Remark 3.1. 

Prom the Proposition 3.1, we know that the function [0,T] 3 t i — > {uj{t))j^'2z 
which belongs to £°°{7L) is derivable ; in particular, the sequence {uj{t))j^2Z is 
bounded uniformly in space and time : 

(3.11) 3M > 0, Vt € [0, T], Vj € I Uj{t) \ < M . 



• We focus now on the total variation of {uj{t)) {c.f. (1.13)). We suppose 
moreover the following restrictions: 
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Property 3.1. 

We fix a real a such that 1 < a < 2. We suppose that the hmiter if related to the 
interpolated values w^, i according to (2.29) satisfies the following inequalities : 

(3.12) ip{r) < (a-2)r, r < 

(3.13) ip{r) < a, r > 0. 



Remark 3.2. 

If (f satisfies the conditions (2. 30). (2. 32), (3. 12). (3. 13) for the particular value 
a = 2 and the complementary restriction : 

(3.14) ip{r) <2r, 0<r<l, 

then the restrictions that we obtain for (p correspond exactly to the "second 
order TVD region" derived previously by Sweby[1984]. For the other values of a, 
(1 < Q! < 2) the estimates (3. 12). (3. 13) are new. We can transcribe the Property 
3.1 and the inequality (3.14) in terms of the ■0-limiter (c./. (2.26)) : 

(3.15) (a) ^(A) > a-2, A<0 
(3.15)(b) ^{J{X) < aX, < A < 1 

(3.15) (c) ^{X) < 2, A > 1. 

Theorem 3.1. Decrease of the total variation 

Let {uj{t)) {j G t e [0, r]) be defined by the method of lines (the hypotheses 
are the ones given at Proposition 3.1). We suppose also that the limiter function 
(f that defines the interpolated values at mid-points according to (2.29), satisfies 
(2.30). (2.32) and Property 3.1 (relations (3.12). (3.13)). Moreover the flux function 
$ of relation (3.2) is supposed to be monotone 

(3.16) ^u{u,v)>0, ^v{u,v)<0 

and absolutely continuous relatively to each variable : 

J Uq J Vq 

Then, the total variation TV(t) of {uj{t)) defined in (1.13) is a decreasing function 
of time : 

(3.17) TV(t + 6') < TV(t) <t <t + d <T . 
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(p(r) 



i 1 




Figure 3.1. Admissible regions for the limiterc/?. The dark shaded region cor- 
responds to the conditions (3. 12). (3. 14) and the clear shaded set completes it to 
describe (2.30) (2.31). 



Proof of Theorem 3.1. 

• We fix some notations : 

(3.18) ipj = (p{rj), Ipj = (p(—), jeTL 

(3.19) ^i+i = ^j+i ~ 5 J £ ^ • 
Then (2.29) admits the following form : 

(3.20) (a) w7^i = Uj + ^(pj A^.+ i 

(3.20) (b) w+^i = Wj+i - ^^j+i ^j+i- 

We also denote by cnj-^ (respectively ) the mean value of (respectively 
over the domain [w7 i , w7, i Ix |w+ i} (resp. {wT, i} x[w+ ^ , w^^, ^ 1) : 

(3.21) $(wT 1 , w+i ) - $(w"_i , w+i ) = a (w" 1 - ) 

J^2 -I 2 J 2 J 2 ■'2 J^2 J 2 

(3.22) $(«7^ , , ) - #(«7^, , , ) ^ p,^, (u^^, - . 

According to the monotonicity (3.16) of the fiux we have : 

(3.23) a^+i > 0, < 0, jeH. 
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When we 



derive formally the total variation (1.13) with respect 



have 

— TV(t) = sgn(uj+i - Uj) -^{uj+i - Uj) 



to time, we 



Then we consider the following sequence 

(3.24) pj^i = -sgn {uj+i - Uj) ^{uj+i - Uj) , j € ^ . 
From the equalities (3.1), (3.2), (3. 20). (3. 22), we have the elementary calculus 



— q; 



/I In / 1 /, 



Then 



A 



and 



(3.25) 



1~ 1~ \ A 



1~ 1~ \ A 



h 



1 - ^( 



1 f 1 - -^j+i + -^j Tj) - a.i Tj f 1 + - 



(l - 5^5^+1 + \vj Ti) - a^_. r, (l + i ^ - i^i+i) _ 
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We integrate (3.25) relatively to time {0<t<T<t-\-9<T), we sum the result 
for j = fc, + 1, . . . , ^ and we factorize 1 and f^j^i. (Abel transformation). 
We obtain : 



(3.26) 



j=k 
j=k 



1 ^3 + 1 



rt+e 1 1 



+ E 



1 + :t T.'Pk-i 1 dr 

2 rfc 2 



1 ) ( 1 - ^V^i+i + ^V'iO- ) 



- / I Afc+ 1 (r) I 3 (^1 - -ifk+i + ^^kVk ) dr . 



The first, second, fourth and sixth terms of the right hand side of (3.26) are 
positive, due to (3.23) and the following inequalities : 

1 ifj+i _ 1 
2 



(3.27) 
(3.28) 



2 rj+i 



• We establish (3.27). We distinguish two cases, according to the sign of rj^i. 
If fj+i > 0, we have, due to (3.13) and (2.30). (2.31) : 

11 1 (fij-Ll 

-^j < -a < 1 < 1 + 

2 2 2 Tj+i 

which proves (3.27) in this case. On the other hand, ifrj+i < 0, we have, due to 
(3.12) : 



■(£>n — -— — < -a 

^^3 n ^ - 2 



5(«-2) 



2^' 2r^+i 

and (3.27) is established. Inequality (3.28) is equivalent , according to (3.18), to 
1 / , 1 1 / 1 \ , 1 



Nonlinear Interpolation and Total Variation Diminishing Schemes 



which corresponds exactly to (3.27). 

• The third and fifth terms of the right hand side of (3.26) are bounded by 

(3.29) 7 = 2M{2 + K)e 

according to (3.11), (3.7), (2.30), (2.31), (2.39) and the constant 7 is independent 
of k and i.. When k tends to —00 and i to +00, the previous inequahties give 
the following estimate : 

(3.30) TV(t) - TV(t + 9) > -AM{2 + K)e. 

By taking t = and 6 <T, this inequality proves that TV(^) is well defined for 
all 9 <T. Moreover, the third and fifth terms of the right hand side of (3.26) tend 
to zero from the previous considerations and the Lebesgue dominated convergence 
Theorem. Taking k — > —00 and £ — > +00 in (3.26), we finally deduce : 

TV{t) - TV{t + 9) > 0, 9>0 

which is exactly (3.17). □ 
Remark 3.3. 

Since the work of Harten [1983], the TVD property is derived from an incremen- 
tal form of the scheme (3.1) Moreover, the inequalities (3.27). (3.28) have been 
previously considered among others by Sweby [1984], Chakravarthy-Osher [1983] 
and Spekreijse [1987]. Here these inequalities are obtained essentially as a con- 
sequence of other considerations concerning the limiter, and not as necessary 
conditions to assume the TVD property. In particular, the limiter function (p{r) 
is not necessarily null for the negative values of the argument r as it is usually 
supposed (see Figure 3.1). To be complete, the sufficient conditions proposed by 
Spekreijse [1987] (see also Cahouet-Coquel [1989]) allow possible non-null values 
for cp{r) when r < but their conditions are not a priori compatible with the 
hypothesis of monotonicity. 



We turn now to an important result of this paper. 



Theorem 3.2. stability. 

Let Uj{t) be the solution of the differential system (3.1). (3. 3); we assume that all 
the hypotheses of the Theorem 3.1 are satisfied. Then we have 



(3.31) inf < Uj{t) < sup uj, j e , < t < T 
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Remark 3.4. 

An immediate corollary of the Theorem 3.2 is the fact that the solution Uj{t) is 
(for each h fixed) defined globally in time, i.e. the time T > introduced at the 
Proposition 3.1 can be taken equal to -\-oo. 

• The proof of Theorem 3.2 contains essentially the ideas given previously 
by Sanders [1983] and Osher [1985]. The new difficulty here is to work with an 
infinite set of indices j (Osher supposed the periodicity in space). So we first prove 
the following property : 

Proposition 3.2. 

For each time t, < t < T, the sequence {uj{t))j^2Z defined by (3.1). (3. 3) has 
a limit M{t) (respectively m{t)) for j — y +oo (respectively j — y — oo) and we 
have 

(3.32) M{t) = M{0), m{t) = m{0), 0<t<T. 
Proof of Proposition 3.2. 

• The existence of M{t) and m{t) is a direct consequence of the convergence 
of the series that defines the total variation for each time t (c./. Theorem 3.1). 

• We now consider the mean value of uj {t) : 

1 ^ 

The function vn derivable and we have 

(3.33) '^ + l(f,,,-f,). 

From the Remark 3.1 and the continuity of the flux relatively to {uj{t)} (see the 
proof of the proposition 3.1), there exists some constant k such that 

Then ^ tends uniformly to zero as TV tends to +oo. Moreover from the Cesaro 
convergence Theorem, v^v tends to M{t) as TV — > oo. We deduce from these two 
facts that M{t) is derivable, and 

(3.35) ^.0. 

The first equality of (3.32) is established. The proof is analogous for m{t). □ 
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Proof of Theorem 3.2. 

• It is sufficient to prove (3.31) for t small enough to establish the property. 
For t G [0,r], we define 

Lbit) = inf u-j(t) . 

Then i — > ii{t) is continuous and we distinguish two cases : 
(i) fi < m. Then for some jo{t), we have 

(3.36) /i(t) = UjM. 

With the notations introduced in the proof of Theorem 3.1, we have : 

du-i ^ , ^ 1 ~ 1 



(3.37) 



/lifjl \ 



and (3.37) is non-negative for j — Jq. Consequently for t' >t and t'—t sufficiently 
small, we have 

(3.38) ii{t') > fi{t), yt'>t, - 1 small enough. 

(ii) ji > m. If the alternative (i) is realized at t = 0, the property (3.38) 

corresponds exactly to the desired result. Otherwise, we have 

(3.39) u(0) = m = inf m° 

and if (3.38) is in default, we consider such that ii{9) < m and r defined by 

(3.40) r = sup {t, t < 6* , ^JL{t) > m } . 

We have < r < ^ and //(r) = m by continuity. According to (3.38), is 
non-decreasing inside the interval [r, 6\ . So that 

(3.41) ii{e) > fi{r) 

which establishes the contradiction. The proof concerning the maximum bound is 
analogous. □ 



4. Convergence, Order of Accuracy 

In this section the parameter h > tends to zero. Two questions are naturally 
related to this problem : 1) What is the order of the truncation error of the scheme 
(3.1) (3.2), (2. 27). (2. 29) as h tends to zero ? 2) Does the function Uh{x, t) defined 
by 

(4.1) Uh{x, t) = u^{t) , {3-\)h <x <{3 + \)h 
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(where the superscript h recalls the dependence towards h of the solution uj (t) of 
the method of lines) converge to the unique entropy solution of (1.1). (1.2) ? The 
following property answers the first point. 

Proposition 4.1. 

We suppose that the entropy solution u{x,t) of (1.1). (1.2) is regular around some 
point (xq, to). Let Vj'{t) be defined by 

(4.2) v'^{t) = u{jh,t), jeTL, h>0, t>0. 

The residual of the semi-discrete scheme (3.1) (3.2), (2. 27). (2. 29), i.e. the quantity 
obtained when we apply the scheme (3.1) to the exact "nodal" values Vj'{t) : 

(4.3) plit) = ^ - rf^K^i'<0 - 1'^^ 

is second order accurate in space : 

(4.4) p';{t) = 0(/i2), /i^O, 3jo, xo = joh 
if the limiter function Lp of the relation (2.29) satisfies : 

(4.5) (a) (p{l) = 1, (p derivable on both sides of +1 

(4.5) (b) (^(-l) + (^(3) = 2 

and if the monotone flux function <l> is absolutely continuous as in Theorem 3.1. 
Proof of Proposition 4.1. 

• The interpolated values i are constructed from the Vj 's according to 
the relations (2.27). (2.29), i.e. 

(4.6) rj , _ r,,h " A . , 

(4.7) (a) v-^,{t) = v] + i(^(r,)A^.+ i = + A_,.+ i 



(4-7)(b) .;_.(t) = - ^^(^) A^_. ^ - 1^, A^._. . 

We develop p^ii) around (a:o, to) assuming that is a mesh point (see (4.4)) 
(4.8) xq = j h for some j E 7L . 

Two cases can occur, according to the nullity of the gradient ^{xq, to) . 
du 

• = —{xq, to) ^ 0. We develop the numerical flux 



dx 



(4.9) /j.+ i = ^{vj + ^ A^.+ i ,Vj + (l-^ ^j+i) A^.+ i ) 
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by the Taylor formula, with the notations 

Oij = = Vj,V = Vj) , = — (U = Vj,V = Vj) , fj = ^{Vj , Vj) . 

We get 

(4.10) ! = + aj^jA.^. + I3j (1 - A,.+ . + 0(A2^.) 

We must now specify the behavior of the (p^s as h tends to zero. We have on one 
hand 

1 u'' 

(4.11) rj = = l-h-^ + 0{h') 

2 u'. 



u 



with = -^^{xq, to) and on the other hand 



1 u" 
(4.12) — = l + /i-^ + 0(h'^) 
rj u- 



When we substitute these developments into (4.10), with the hypotheses (c./. 
(4.5)) 

(4.13)(a) ^(1 + = 1 + + O(e'), e^O, e>0 
(4.13)(b) ^(1-0 = l-^^- + 0(a, ^^0, ^>0. 



We obtain (assuming ^ > ; the other case is analogous) 

3 

\. . , . . 1 



(4.14) 



/j+i = U + ^^^(«j + /5j) + ^^^^'(-«i<^- - + 

+ ^^;'(c., + /5,) + ^(^)V; + 0(/z^) 



(4.15) = {^uu^^uv^^vu^^uu){Vj,Vj) = -j:i^{Vj) 



where /j' is defined by the formula : 
Similarly, we easily obtain 

.,h _ „,h // 

(4.16) r,_i ^ = 1-^^ + 0(/i2) 

u'- 

(4.17) (a) v'j-i = l-/i^^'_ + 0(/i2) 
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u'' 



(4.17)(b) = l + + 0(/i') 



We insert the latter expressions in the fohowing expression of the numerical flux 
at {j - : 

(4-18) = ^{^j + {l^j-i - l) , Vj - A._i ) 

and we obtain finally: the other case is analogous) : 

(4.19) 



+ ^<(a,+/3,) + ^(n;f/; + 0(/z^). 



(4.20) (a,+/3,)n;. = f\vj{t))u'^ = -^^Uj{t) 



We subtract (4.19) from (4.14) and we also notice that 

d_ 
di 

Finally we get 

(4.21) p^;{t) = 0{h') 

as desired. 

• We consider now the case of an extremum : u', = -r— (iCn , to) = 0. Then 
— (a:o , to) = due to the conservation law (1.1). We focus on the development 
of the ratio rj of the slopes : 

In a similar manner, we have: 



1 _ fwy + o{h) 

r~^i ~ \u'; + 0{h) 



(4.23) — = I ], , = 3 + Oih) . 



Then the flux at the interface (j + |) admits the following form : 

(4.24) . = + \h^v!l (ia,- v^(-l) + /3,- (1 - i^(3))) + Oih") . 

The expressions r^-i and ^ are computed according to (4.23) and (4.22) respec- 
tively. Then the numerical flux at (j — |) /i can be developed as 

(4.25) = - \h^v!'^ ((1 - \p{Z))a^ - \p{-\)^i) + 0{h^) 
and we finally get in this case : 
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that is 

(4.26) p^(t) = -\hu'! (^ftll±£B _ i) (a. _ + o(ft^) 

and the proposition is estabhshed. □ 

Remark 3.1. 

In a previous paper, Osher-Chakravarthy [1984] claimed that if a semi- discrete 
scheme which takes the form 

with 

(4.28) i > 0, i > 0, jeTL 

is necessarily first order near a non sonic critical point (|| = 0, f'{u) ^ 0). The 
numerical semi-discrete scheme (3.1)(3.2), (2.27). (2,29) studied herein admits the 
form (4.27) (4.28). More precisely, we have, according to (3.20). (3.22) : 

1 1 

(4.29) (a) C,+ 1 = -^^.+ 1 ( 1 - -^,+1 + -^,r, ) 



(4.29) (b) 




Near a critical point, the expressions (3.27) and (3.28) vanish with h and according 
to (4.5), the factors Aj_^i are 0(/i^). Moreover the coefficients and ctj_i 

are bounded. The second order is clear. In fact the consistency property assumed 
in the previous reference (formula (2.11)) is not relevant for a critical point while 
the term in (2.10) is null. The result can also be stated intuitively : near 
a regular extremum, the curve looks like a parabola and its interpolation by a 
constant is second order accurate in that case. 

Our condition (4.5) (b) which is the key for this second-order accuracy has been 
independently proposed by Wu [1989]. 

• We study now the convergence of the sequence Uh{x, t) defined in (4.1) 
to the entropy solution of (1.1) (1.2). We have not proved that the scheme studied 
previously converges to the entropy solution ; therefore we give known results 
concerning the convergence of Uh{x, t). First we have the following property : 

Proposition 4.2. 

The sequence of functions Uh{», •) defined for (x^t) G IRx]0, -|-oo[ from the line 
values Uj{t) according to the relation (4.1) are bounded in each compact set in 
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L°° n BV(IR X 0, +00 [) where the space BV is defined according to Volpert [1967]. 
This bound is independent of h > 0. 

Proof of Proposition 4.2. 

In Part 3, we have proved the BV property relatively to the space variable x 
(Theorem 3.1) and the L°° stability (Theorem 3.2). It is now suflacient to prove 
the BV stability for the time variable. Let T and At > be two fixed times. We 
have : 

J ( y dx \ Uh{x, t -\- At) — Uh{x, t) I ^ < 

pT /.At o 

< hJ^T— de(\ Uj+i - I + I Uj - \){t + e) 

where the constant C has been introduced in (3.7). Therefore the BV property in 
time follows from the BV stability in space (Theorem 3.1). □ 

• From Proposition 4.2 and the compactness of the injection L°° fl BV — > 
h°° (see Volpert [1967]) we deduce that a subsequence of uh converges to a weak 
solution of (1.1) (1.2) but not necessarily to the entropy solution. To prove that 
a subsequence of uh converges to the entropy solution it is suflacient to estab- 
lish a discrete entropy inequality (Harten-Hyman-Lax [1976]). Recently, Osher 
[1984] has proposed a sufficient discrete entropy condition to establish the discrete 
entropy inequality : 

(4.30) / r]'{w) ( /^.+ 1 - fiw) )dw < 0, j e 7Z . 

When /j+i is computed according to the MUSCL method {e.g. the relations 
(3.2), (2.27)-(2.29)), Osher [1985] proved that when the flux function $ in (3.2) 
is monotone (which has been always assumed in this paper), the sequence Uh 
converges almost everywhere to the unique entropy solution of (1.1) (1.2) if the 
limiter function Lp{m) satisfies the monotonicity property 2.4 and the following 
severe restrictions : the function (^(•) is null for negative arguments and the 
function (/p(») satisfies the following property : 

(Ao^\ J < 2 max 0,min( I 

l^-'J-L j \ L ^Uj — Uj+i Uj — Uj + i / J 

if Uj > Uj+i , for each j G 2L 
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with the following definition for the mean value at {j-\-^)h 
(4.32) f(u.^.) = I f{w)dw 



The condition (4.31) couples the interpolation at the interface values and the 
(physical) flux function /. This result is at our knowledge the best one concerning 
the convergence of the sequence Uh towards the entropy solution of the conservation 
law if one uses the MUSCL approach. For other studies concerning this problem 
we refer also to Osher-Tadmor [1988] and Vila [1988]. 



5. Time Discretization 

We focus now on the discretization in time of the semi-discrete scheme (3.1). 
(3.2) (associated with (2. 27). (2. 29)). We assume that this semi-discrete scheme 
is Total Variation Diminishing, i.e. that the sufficient conditions (2.30)- (2.32) 
and (3. 12). (3. 13) are satisfied. We consider the following three schemes in time 
(At = t'^~^^ — is the time step and we denote by A the ratio ■^) : 

First Order Scheme 

(5.1) - «») + ^(*(«;^ , - , ) = o. 

Modified Euler Scheme (Heun scheme) 

(5.2) (a) ^^(ur - u]) + , - , ) = o 

(5.2) (b) ' - uj) + £ ( (/jv, + - (/;_, + f;^}) ) = 

with the fluxes given by the relations 

(5.3) f^, = , , 1;:^ = - . a;;:- . 

Predictor- Corrector Scheme 

(5.4) (a)) A _ „») + J- ( . - z;^ . ) = 

(5.4) (b)) - nj) + 1^ ( f;^^ -Cl) =Q. 
As in the previous cases we have the same type of notations : 

(5.5) /•;. = , , d = $(«;:r , . 

In the following, we study the TVD property of the scheme (5.1). (5. 4) under some 
ad hoc restrictions. 
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Property 5.1 

We suppose that the hmiter function defined in (2.9) satisfies the following 
inequality : 

(5.6) 3M>0, Vr>0, ^ < M. 

r 

Prom the Figure 3.1, we have clearly M > 1. 
Proposition 5.1 

The first order scheme (5.1) joined with (2. 27). (2. 29) is TVD under the hypothe- 
ses (2.30). (2.32), (3. 12). (3. 13) and (5.6) for the limiter function if the following 
Courant-Friedrichs-Lewy condition holds : 

(5.7) ^(|a»^,l + |/3;^,l) < jeTL 

with the partial gradients of the fluxes otjj^i and /^j+i deflned by the formulae 
(3.21). (3.22). 



Proof of Proposition 5.1 

• The incremental form of the scheme (5.1) has the expression 



1 1 



1 



3 

(we drop the index n when there is no ambiguity) . Therefore it is under the form 
considered by Harten [1983] : 

(5.9) = + C^.+ . A^+. - A^._. 

with 

(5.10) (a) c^+i = A (-/?^+ . ) ( 1 - hfi+i + Im ) 

(5.10) (b) D,_. = Aa,_. (l + i^-i^i-i). 

Harten's conditions must be satisfled : 

(5.11) (a) C,.+ i > 0, (b) D.^. > 0, (c) C^.+ i i < 1. 

• The inequalities (5.11)(a)(b) are a direct consequence of the Theorem 3.1. 
On the other hand, we have : 

(5.12) (a) l + < 1 + y 
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(5.12) (b) l--^.+i + -^.r,- < 1 + Y 

and the inequality (5.11)(c) is a direct consequence of the CFL condition (5.7). □ 
Proposition 5.2 (Shu-Osher [1988]) 

The modified Euler scheme (5.2). (5.3) is total variation diminishing under the 
same hypotheses as in Proposition 5.1 and the same CFL restriction (5.7). 

• The predictor-corrector scheme is more complicated to study. We begin 
with the following lemma : 

Proposition 5.3 

We consider a numerical scheme given by the following incremental form : 

(5.13) u]+^ = Uj + a^.+ 3 A^.+ 3 + c,+ i A^.+ i - dj_i A^._ i - h^_^ Aj_3 . 
This scheme is TVD if the following inequalities hold for each integer j : 



(5.14)(a) 




(5.14)(b) 




(5.14)(c) 


1 ^n-i- ^ d A _i_ 1 
J 1 2 J ^ 2 


(5.14)(d) 


d^x 1 — e„_i_ 1 > 

2 2 


(5.14)(e) 





> 



This proposition generalizes the one proposed by Harten [1983] and used previously 
in this paper. The proof is straightforward. We omit it. 

• We now analyze the TVD property for the predictor-corrector scheme (5.4). 
(5.5). This scheme does not belong to the general class considered by Shu-Osher 
[1988]. In fact, we have not proved that under the form (5. 4). (5. 5) the predictor- 
corrector scheme is total variation diminishing but no conter-examples have been 
found. Nevertheless we modify the scheme (5.4). (5.5) in the following way : 

Definition 5.1 

For k equal to n and n-\- ^ let Cj_^i and Dj_i be defined by the relations 

k \ V J+2 J+2 / V J+2 3-2 



(5.15)(a) C'. = -A 

2 



(5.15)(b) D''_^ = A ^ ^-^-^ 

3 2 Uj+i-Uj 
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We define two versions of the modified predictor- corrector scheme by the incre- 
mental relations 

(5.16)(a) =u, + \{c^^, A,+ . - A,_ . ) 

(5.16) (b) = + 1 ( . a;;| - D^_^ a;:| ) . 

If k is equal to n (respectively n + |) in relation (5.16) the predictor-corrector 
scheme is explicit in time for the two steps (respectively implicit for the predictor 
step and explicit for the corrector). Notice that in equation (5.16) the index k is 
the same in both steps (5. 16) (a) and (5.15)(b) whereas in the original predictor- 
corrector scheme (5.4) (5.5) k is equal to n for the predictor and k is equal to n + 1 
in the corrector. We have the following property : 

Proposition 5.4 

The predictor-corrector schemes defined in Definition 5.1 are TVD under the hy- 
potheses of Proposition 5.1. 

Proof of Proposition 5.4 

• We first develop and we drop the exponent n for all the variables 

defined at time nAt and the exponent k of (5.16) associated with the coefficients 
C and L> : 

(5.17) A;;| = (l-i(C,+ . A,+ . + ic,+|A,+| + li?,_.A,_.. 



• We substitute (5.17) in equation (5.16)(b). We obtain in this manner : 




This representation is of the type proposed in (5.13). The inequalities (5.14) take 
the form : 



(5.18)(a) C,.+j Cj_. > 

(5.18)(b) C,+ . (1 - i(C,+ . + D,^,)) - + C,_.)C,+ . > 
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(5.18)(c) 



1 



(5.18)(d) (C^^.+D^^.) (l-^(C,+ .+D,+ .)) - \d^_, (C,+ , > 
(5.18)(e) > 0. 

From the inequahties (5.11) the inequaUties (5.18) (a) and (e) hold clearly. On 
the other hand, we have : 

left hand side of (5.18)(b) > C^+i (l - ^) - ^C'j+i = 

and (5.18)(b) is established. The inequality (5.18)(d) holds by the same argument. 
Finally, the left hand side of (5.18) (c) is dominated by (C'^+i + -C^j+|) and this 
last inequality is established according to (5.11)(c). □ 



Remark 5.1 

As mentioned previously, we have not proved that under the natural form (5. 4). (5. 5) 
the predictor-corrector scheme is TVD. In fact, the inequalities of the type (5.18) 
contain in that case a mixing of the Cs and i^s with time steps n and n + |. The 
simple algebra derived in the last proposition is not possible anymore. We also re- 
mark that Definition 5.1 proposes a scheme under its incremental form. Therefore 
this scheme is not a priori conservative. These two remarks illustrate the difficulty 
in maintaining both the conservative form of the numerical scheme and the TVD 
property for a predictor-corrector scheme constructed from the first order scheme 
(5.1). 



• In the last part of this section we study the order of accuracy of the predictor- 
corrector scheme (5. 4). (5. 5) and the modified Euler scheme (5. 2). (5. 3) in the sense 
of the truncation error near a regular local extremum of the solution (t, x) i — )■ 
u{t, x) of the conservation law (1.1). We recall that in Proposition 4.1 we have 
proved that the truncation error (4.4) for the method of lines is second order 
accurate in space if the limiter function (f satisfies the condition 

(5.19) ^(-l) + (/p(3) = 2. 

We first make the following two hypotheses concerning the numerical flux function 
$ of relation (3.2) and the limiter function <^(— 1) for the interface interpolations. 

Property 5.2 

We suppose that the monotone flux function $ = $(w, v) satisfies the condition : 
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(5.20)(a) df{w) > ^ — (w;, w) = 

(5.20) (b) d/H < ^ -^iw,w) = 0. 

• We note that the monotone flux functions proposed by Godunov [1959] 
and Engquist-Osher [1980] satisfy the hypothesis 5.1 whereas the numerical flux 
associated with the Lax-Wendroff [1960] scheme : 

(5.21) ^^{u,v) = ^(f{u) + f,u) + ^(f{v)-^iv) 

(where fi is chosen such that d/(w;) < — fi in the domain of variation of w) does 
not. 

Property 5.3 

The limiter function r i — > (p{r) satisfies the identity 

(5.22) (1 - a) ^(^^) + (1 + ^) ^{^) ^ 2 , Va € [0, 6] 

for some parameter 6 > and the discrete relations 

(5.23) (A;+l) ^(^) - (A;-l) ^(^) = 1, /c = -2, -1,0, +1, +2. 

Proposition 5.5. Accuracy of the predictor-corrector scheme 

If on one hand the flux function $ is regular {i.e. absolutely continuous relatively 
to each variable, monotone (relations (3.17)) and satisfies the property 5.2, and 
if on the other hand the limiter function if satisfies the monotonicity property 
2.4, the property 5.3, and is derivable at the particular value r = 1, then the 
predictor- corrector scheme (5.4). (5.5) admits a truncation error 

(5.24) pih, u) ^ ^(^(^,-, r+i) -n(a;„ T)) + ^ (/;+|(h) - f^^Hu) 

of second order accuracy : p(/i, u) = 0(/i^) if the following CFL restriction is 
satisfied : 

At 

(5.25) ^ sup \df{w) \< 5. 
Proof of Proposition 5.5. 

Three cases are distinguished by the fact that the first derivative of u is null or 
not. We use the same notations as in Proposition 4.1. 

= ^{xj, t") ^0. We have 
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(5.26) 



= 1- 



u 



n _ n 
^j+k ^j+k-1 



XX 



h + 0(/i^) , j, k integers. 



When we introduce this relation into (4.7) we have 

n, lb 



(5.27) 

and 
(5.28) 



V 



j+k+. 



- i At /, - i At h(k+^) + 0{h^) 



V 



'j+k 



1 



Vj -\- khux -k'^h^ Uxx — Ti^t fx — -Atkh fxx + 0(^'^) 



1 



1 



Tl-\- — 1 

So that the ratio rj_^^ at t = (n+ ^)At also satisfies the relation (5.26) (!), i.e. : 



(5.29) 



j+k 



We insert again this relation into (4.7). We deduce 



h + 0{h^) 



(5.30) 



V 



n, ± 
3+k+ 



1 = Vj + {k+-)hux + -(k^ + k+--ip'{l))h^ 
-^Atfx - ^Ath(k+^) fxx + 0(/i^) 



and the numerical flux can finally be developed at the point {j + k + as 



3+k+,^ 
''j+k+h 



(5.31) 



= fj + (k+l)hfx - lAtf'fx + (k+^)'h'f 
- i At (/c + ^) /i^ (/' /^^ + /' ux) e f, 

+ iAt2 [f\f"Uxfx + nxx) + f'ifxf] + Oih') 



XX 



XX 



Then the finite difference of fluxes at the point Xj is given by : 

(5.32) /;+^ - /;+^ = fxh-l-Ath (f fxx + /" Ux fx) + 0{h^) 
Finally the truncation error is evaluated by 
(5.33) 



p{h, u) = ut + ^Atutt + fx - I {f fxx + f Ux fx) + 0{h^) 



and the result is in this particular case a direct consequence of elementary algebra 
from the partial differential equation (1.1). 

• Ux = , Uxx 7^ . We observe first that v"^^^ is defined by the interme- 
diate values for k = —2, —1, 0, +1, +2 because the interpolation (4.7) is 
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a three-point scheme. Moreover these nodal values at t = (n + |)At are defined 
from the given values Vj'_^f^ = f^) for k = —3, —2, —1, 0, +1, +2, +3. We 

have : 



(5.34) 
(5.35) (a) 
(5.35) (b) 



j+k 



k 2 



f + 0{h) 



V 



j+k+i 



Vj + ^h^Uxx 



n, + _ 



W: ; , 1 = v.- + - /l iX^a; 



(fc+lf -(fc+i)v(^)]+0(ft^). 



Therefore the flux difference (^/f, l. = i, , i — ^ i are evaluated by the de- 
velopments 



(5.36) 



f7+k = kh^f u,, + 0{h^) , /c = -2, -1, 0, +1, +2 



because the coefficients of the partial derivatives of the flux function $ in Taylor 
formula are null according to the property 5.3 (relations (5.23)). The values at 
the intermediate time admit the development 



(5.37) 
and we have 
(5.38) 



v:t^ = Vj + ^k^h^u,, - ^Atkh'^fu,, + 0{h^) 



r 



n+i 



-I + 2k - Xf 



1 + 2/c-A/' 
Using the Taylor formula again, we finally get 



+ 0{h^), /c = -l, 0,+l 



(5.39) 



2 J 1 



2 



du 



3 + A/' 

+ A/ 

3-Af 

l-Xf 

+ 2^-'-^^^^ItTaF 



+ A/'. 



+ 



+ 



To fix the ideas, we suppose /' > 0. Then ^ is null according to Property 5.2 
of the numerical flux function. Moreover, the coeflacient in front of ^ is also 
null thanks to Property 5.3 (relation (5.22)) if we have Xf'<5 which is exactly 
the CFL condition (5.25). The proof is analogous when /' < 0. The sonic case 
/' = is very simple because the brackets of (5.39) are both null according to 
(5.19) (which can be seen as a particular case of both (5.22) and (5.23) with cr = 
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k = respectively). The truncation error (5.24) is again evaluated according to 
relation (5.33) and the property is established in this second case. 

• Ux = , Uxx = 0. We have clearly 

«;;ti = V, + om 

because the function (/?(•) is bounded and is 0{h^). It is straightforward 

to develop the truncation error (5.24) and to verify the property. This ends the 
proof. □ 

Remark 5.2 

If the limiter function satisfies Property 5.3 and is affine, it takes necessarily the 
form 

(5.40) ip{r) = 1 + a(r- 1). 

This kind of function is not compatible with the monotony- convexity conditions 
(2.30)- (2.32) and the sufficient TVD conditions (3. 12). (3. 13) if (5.40) holds for ah 
the real values of the argument r. Therefore the actually known limiters {e.g. 
minmod, superbee (see Roe [1985]), MUSCL (see Van Leer [1977]), etc. ) do not 
satisfy both (5. 22). (5. 23). Nevertheless , the function defined as follows 



(5.41)(a) 


(p{r) = 


0, 




r < -3 


(5.41)(b) 


(f{r) = 


1 . 


+ 3), 


-3 < r < - 


(5.41)(c) 


if{r) = 


r 

"2 ' 




-1 < r < 


(5.41)(d) 


ip{r) = 


5 r 

T' 




0<r<l 


(5.41)(e) 


cp{r) = 


1 . 


+ 3), 


i<r<3 


(5.41)(f) 


(f{r) = 


3 

' 




r > 3 



is afRne and satisfies (5.40) for r e [—3, —1] U [|, 3]. Then (5.23) is satisfied. 
Moreover (5.22) holds for 6 = ^ and Property (5.1) is also satisfied with M = |. 
This limiter satisfies also the TVD conditions (3. 12). (3. 13) of Property 3.1 with 
q; = I and the properties 2.4 and 2.5 of monotonocity and convexity. Finally we 
observe that the interpolate value 

= ^3 + \ '^(^i) (^i+i ~ ^j) = L{uj-i , Uj , Uj+i) 

(according to the notation proposed in (2.1)) corresponds to the Lagrange polyno- 
mial interpolation of degree 2 if and only if if satisfies (5.40) with a = \. In the 
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following we will denote this new limiter function (5.41) by "Lagrange limiter". 
We note also that the idea of using Lagrange interpolation for the construction 
of limiter functions was at our knowledge first proposed by Leonard [1979] with 
the so-called Quick limiter. For a modification of this limiter that does not sat- 
isfy the monotonicity condition for negative values of the argument r, we refer to 
Cahouet-Coquel [1989]. 



(p(r) 




Figure 5.1. The Lagrange limiter (defined by the relations (5.41)). 



Proposition 5.6. Accuracy of the Heun scheme 

We assume that the flux function $ and the limiter function satisfy the same 
hypotheses as in Proposition 5.5. Then the truncation error 



1 



(5.42) 



+ 



1 



2 Ax 



n+l , 



'3 + 



'3 + 



of the modified Euler scheme (5. 2). (5. 3) is of second order accuracy if the 
following more restrictive CFL condition is satisfied : 



(5.43) 



At 

Ax uu 



f'M\< 



Proof of Proposition 5.6 

As in Proposition 5.5 we divide the proof in three steps, following the value of 
and Uxx- 
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• Uxixj , t"^) 7^ 0. We develop the numerical flux at time and x = {k -\- 
j -\- as in Proposition 5.5 : 



(5.44) 



On one hand, we have 

(5.45) 5fJ = hf, + 0{h^) 
and on the other hand : 

(5.46) u]+l = Uj + {khu:,- l^tU) + ]^h{ehu^^-2k/^tU^)+0{h^). 

Therefore the ratio of gradients after the "tilda" step (5.2) (a) is given by (5.26) or 
(5.29). The extrapolated values at the interface {j + k+^) are again evaluated 
according to (5.30) except that At must be replaced by 2 At. In an analogous 
way, the tilda flux defined in (5.3) admits a development (5.31) with At. replaced 
by 2 At. We deduce : 

(5.47) Tf^ = hf,- Ath(f' f^^ + f'u^f,) + 0{h^). 

The truncation error admits exactly the form (5.33) and the property is proved in 
this case. 

• Ux{xj, V^) = , Uxx{xj, t") 7^ 0. The diff^erence of fluxes around the point 
{j + k)h is given according to relation (5.36) due to Property 5.3. We have : 

(5.48) r,,. _ 1 ^ _ 2,;, + Oih) 

according to (5.43) and (5.22). Therefore the development of the flux difference 
at point j h and at time (n + l)At is simply developed : 

(5.49) S/f ' ^ r+l - = /7xx + 0(/z3) 

and the truncation error admits flnally the form (5.33). □ 



6. Numerical Tests with the Advection Equation 

In this section we present numerical comparisons of the schemes studied in Part 
5 for a simple model problem. We show that the truncation errors studied at 
propositions 5.5 and 5.6 are numerically second order accurate for the Lagrange 
limiter (5.41) especially at a nonsonic extremum. We study also the effect of the 
restrictive CFL condition (5.43) on the accuracy of the two step Heun scheme. 



FRANgois Dubois 



Moreover we compare the accuracy measured in the discrete norms 
and in terms of the truncation error (see also Harten et al [1987]). 

We focus on the one-dimensional linear advection equation 
Ou Ou 

(6.1) — +a— =0, 0<a:<l, t>0, a = l 

at ox 

with periodic boundary conditions : 

(6.2) w(0) = w(l). 

We look at the advection of a regular profile with a (local) extremum at x = ^ 
without annulation of the third derivative at this point (which is the case with a 
sinus-type profile) giving unexpected superconvergence for the truncation error at 

X — 2 '■ 

We define our profile as a class periodic function on the interval [0, 1] 
according to the conditions 

(6.4)(a) = 1 - (x - i)' + (x - i)' , \<^<l 

(6.4)(b) uyo) = '^iO) ^ 2^(0) ^ .»(!) ^ '^(1) ^ 2^(1) = 

1 3 

(6.4) (c) polynomial of degree 5 in [O, -] and in [-, l] . 

The uniqueness ofw° defined by (6.4) is clear due to the use of the classical finite 
element interpolation in the intervals [0, |] and [|, 1]. 

• We study the Cauchy problem (6.1) (6.2) associated with the initial condition 

(6.5) w(0, x) = u^{x) , < a: < 1 

during one period, i.e. for time t in [0, 1]. At time t = 1, the exact solution 
of (6.1). (6. 2). (6. 5) is equal to the initial condition and the error is therefore 
easy to evaluate. We discretize the equation (6.1) on meshes with cell diameters 
Ax = h that are constant on each mesh and chosen as powers of | : 

(6.6) h = 2-^ /c = 3, 4, 5, 6, 7, 8, 9, 10, 

and the computational nodes Xj (the cell centers in the finite volume approach) 
are located at 

(6.7) Xj = {j-l)h, j = 1, ...,Nih) = i. 

We have tested three types of temporal schemes : the so-called "five points TVD 
scheme", the predictor- corrector scheme, and the Heun scheme (second order 
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(6.io)(b) =ui-o («;+%^-- - 



— I ~n+l,—\ ( — I ~n+l, 



i+1/2 ' ^i+1/2 ; V"j-l/2 ' "'i-1/2 



Runge-Kutta) . All these schemes are parameterized by the interpolation u^_-^i^ 
(at the left side of the interface Xj_ii2) from the cell center values (6.7) in the 
spirit described in Part 2 of this paper and detailed in the following. The temporal 
schemes are associated with the Courant number 

(6.8) a = a^. 

To be precise we give the formulae for this temporal integration. 

• Five points TVD scheme (see e. g. Sweby [1984]) : 

(6.9) = Uj - a a{uj - Uj_i) + {1 - a) {uj_^^^^ - ^7-1/2) 

• Predictor- Corrector Scheme (see also (5. 4). (5. 5) ) : 

(6.10) (a) = Uj - ^ (^7+1/2 ~ V1/2) 

n+l/ 
S-1/2 

• Heun scheme (see also (5. 2). (5. 3)) 

(6.11) (a) u]+^ = Uj - (t(w7^^/2 - V1/2) 

(6.11) (b) u]+' = - f 

The interpolation at point is chosen according to a three point interpolation 

as in the previous sections : 

(6.12) = + 2'^\ uj+i-Uj ) V^J+i"''^') 

or with the UN02 interpolation proposed by Harten-Osher [1987]. Concerning the 
limiter functions we have restricted ourselves to the following ones : 

(i) First-order upwind scheme 

(6.13) (^^P^i'^d(r) = 

(a) hax Wendroff type scheme 

(6.14) (f^^ir) = 1. 

We remark that (6.14) joined with (6.12) defines the Lax- Wendroff [1960] scheme 
only when it is associated with the temporal scheme (6.9). The two other cases 
justify the "type" restriction in our denomination. 

(in) Van Leer's MUSCL limiter [1977] 

(6. 15) (a) v'^uscL^^^ = if r < 0, ip^^^^^{r) = 2rif0<r<i 
(6.15)(b) ^^uscL(^^ ^ .f 1 < ^ < 3, ^MUSCL(^) ^2 if r > 3. 
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(iv) Min-mod limiter (Harten [1983]) 

( ^^''^^°<^{r) = if r < 0, 

(6.16) I ^minmod^^j = T if < f < 1 , 
[ (^minmod^Y) = 1 if r > 1 

(v) Superbee limiter (Roe [1985]) 

(6.17) (a) (^^^P«'^B(r) = if r < , (^^"P^^^(r) = 2r if < r < ^ 
(6.17)(b) ^^^P«'^^(r) = 1 if ^ < r < 1, 

(6.17) (c) (^^^P«^B(r) = r if 1 < r < 2, (^^^P«^B(r) = 2 if r > 2 . 
(^m^ The Lagrange limiter proposed in Part 5 (relations (5.41)). 

(vii) Min Mod absolute value limiter, that we define according to 

(6.18) ip^^^ir) = I ^ I ^ _ 

(via) The UN02 interpolation, which is a five point interpolation that is defined 
according to 

(6.19) minmod(a, 6) = a (^"^^""^"^ (-) 

(6.20) d'^Ujj^xji = niinmod ( Uj-\ — 2uj + i '^j ~ '^j+i + ) 

111 

(6.21) u7^^/2 = + 2 ^^%-i/2 , Uj+i-Uj-- ) . 

We note that this UN02 scheme has a stencil of a priori 7 points for the one step 
temporal scheme (6.9) (only 5 points in the interpolations (i)-(vii)) and 13 points 
for the temporal schemes (6. 10). (6. 11) (compared to 9 with the other schemes 
described herein). 

Known results established in Sweby [1984] are precised by the following proposi- 
tion : 



Proposition 6.1 

The discretization of the scalar advection equation (6.1) joined with the five point 
scheme (6. 9). (6. 12) and some limiter function (p admits a truncation error 



\U(Xj, t 



n+l 



(6.22) p{h,u) = < 



Ax ^ ^ 
h\Ax J \ 



V 



J -1/2 



{u, n 
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(where '^^j_^i/2(^5 is a notation for interpolation (6.12) applied with the val- 
ues of the solution interpolated at the grid points a;j_)_i/2) which is second order 
accurate 

(6.23) p{h, u) = 0{h^) 

at a local extremum of the function u : 

(6.24) —{x = Xj,t = e) = 0, \fh 

if the following condition (proposed in Part 4) holds : 

(6.25) ^(-1) + ^(3) = 2. 



Proof of Proposition 6.1 

• From the relations (5.35) (a) we have 

. - - 1,2 

(6.26) - = 2 ^ ^ 

and p{h, u) admits the development 



(6.27) p{h, u) 



du 1 d'^u 1 2 7 

)2. 



1 f)^ii /I 1 \ 

Thus estimate (6.23) is a consequence of (6. 25), (6. 8) and (6.1). 



□ 



We can compare on figure 6.1 several profiles obtained for different meshes at time 
t = 1 with the Heun scheme. In order to quantify the errors, we set 

N{h) 

(6.28) 11 UH{t = l)-u{t = 1) 11^1 = hJ2 \ ^j,hit = 1) - u'iooj) 



(6.29) \\uhit=l)-uit=l)\\^2 -- 

(6.30) \\uh{t=l)-u{t=l) ll^oo 



Nih) ^ 1/2 

h J2 \^j,hit = l) - u%xj) 

j=l, ...,N{h) 
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advection profil regulier, schema de Heun, limiteur de Lagrange 
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Figure 6.1. Advection of a regular profile. Solutions obtained with the La- 
grange limiter and the 5 points uno2 scheme. 
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• We present four numerical experiments associated respectively to the fol- 
lowing choice concerning the temporal scheme : 

exp 1 Five points TVD scheme (6.9) CFL : cr = | 

exp 2 Predictor-corrector scheme (6.10) CFL : cr = | 

exp 3 Heun scheme (6.11) CFL : cr = | 

exp 4 Heun scheme (6.11) CFL : a = ^. 

For each of these numerical experiments we have advected the initial profile defined 
in (6.4) during one period with seven diff^erent meshes defined in (6.6) and the eight 
different fluxes (i) to (viii). We present on flgures 6.2 to 6.5 the errors 
defined in (6.28)- (6.30) as functions of log2(l//i) and the truncation error (deflned 
respectively in (6.22), (5.24) and (5.42)) at the particular point x = ^ where the 
maximum occurs (it is always a mesh point). All the residuals are proposed in 
terms of their base 2 logarithm. Therefore if some residual r{h) is "at the order 



a" 



(6.31) r{h) = C/i" (l + o(/i)) 

the slope of the curve log2[r(/i)] as a function of log2(l//i) is simply equal to 
—a and is directly readable on the figures. We notice first that the truncation 
errors at the extremum {x = |) have exactly the behavior predicted in (6.23), 
(5.25) and (5.43); with c" = ^ the order of the truncation error at the extremum 
is only 1 with the Heun scheme (Fig. 6.4.d) therefore it is of order 2 with the 
same numerical scheme and c" = ^ (Fig- 6.5.d). The five point scheme (Fig. 
6.2) show that the MUSCL limiter and the UN02 scheme are really second order 
accurate in the L°° norm (Fig. 6.2.c). When we use a two step scheme and a 
Courant number c = | (Fig- 6.3 and 6.4) the Lagrange limiter is "better than 
second order accurate" (order 2.2 in the L°° norm) and the other limiters are 
second order accurate but "only" in the norms LI and L2. With the Heun scheme 
and a severe CFL restriction (Fig. 6.5) the results of the Lagrange limiter are 
degraded despite the fact that the truncation error remains everywhere second 
order accurate ! 
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Figure 6.2. Advection of a regular profile. Errors for the flux corrected trans- 
port type temporal scheme with Courant number = 0.5. 
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Figure 6.3. Advection of a regular profile. Errors for the predictor-corrector 
temporal scheme with Courant number = 0.5. 
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Figure 6.4. Advection of a regular profile. Errors for the Heun temporal 
scheme with Courant number = 0.5. 
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Figure 6.5. Advection of a regular profile. Errors for the Heun temporal 
scheme with Courant number = 0.25. 
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7. Conclusion 

In this contribution, we have proposed a concept of convexity preserving property 
for the approximation of scalar conservation laws with the Van Leer's numerical 
method. New restrictions have been derived for the construction of the associated 
limiter functions. We have proved that the method of lines is well posed for uq 
in £^ n BV and new conditions on the limiter have been given to prove the TVD 
property of the associated solution. We have proved that second order accuracy 
(in the sense of the truncation error) can be maintained, even at a nonsonic ex- 
tremum ; this last property and Total Variation Diminishing are compatible with 
a discretization in time with a two-level Runge-Kutta scheme under new restric- 
tions on the limiter function that have been detailed. An example of such a limiter 
(the so-called Lagrange limiter) has been proposed and implemented in the case 
of the advection equation. Theoretical results on orders of convergence have been 
confirmed numerically. This concept of convexity preserving interpolation could 
be extended in two directions : the treatment of multidimensional scalar problems 
and the approximation of systems of conservation laws. 
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